chapter8- Regression
alfonso
May 21, 2016
关于“回归”
例如一位工程师想找出对桥梁恶化最重要对影响因素（时长，桥梁设计，材料，天气等），来确定这些关系的数学形式。以下几点会有助于这个过程：
- 找出有关的因素，及其相对重要性
- 对于异常值要谨慎处理，可能会有大发现
- 模型可能会有多个
本章的思路： 1. 如何表示模型 2. 展示各种回归方法 3. 判断回归模型是不是有问题 4. 找出特别的观测值 5. 纠正模型
最小二乘回归
R 中的公式
如何用 R 表示 输入与响应的公式？ Y ~ X1 + X2 + ... + Xk 就表示一个线性回归的公式。
| 符号 | 应用 |
|---|---|
~ |
用来分开输入与响应，左边的位响应 |
+ |
间隔开输入变量 |
: |
两个输入变量间的相互作用 ， y ~ x + z + x:z |
* |
所有有可能的相互作用的简写， y ~ x * z * w <=> y ~ x + z + w + x:z + x:w + z:w + x:z:w |
^ |
标记相互作用的上限， y ~ (x + z + w)^2 <=> y ~ x + z + w + x:z + x:w + z:w |
. |
data frame 中除去依赖变量的所有变量的占位符， y ~ . <=> y ~ x + z + w |
- |
从方程中移除某项， y ~ (x + z + w)^2 - x:w <=> y ~ x + z + w + x:z + z:w |
-1 |
去掉截距，例如： y ~ x -1 强制回归线通过原点 |
I() |
表示括号中的需要预先处理，例如：y ~ x + I((z + w)^2) 将会扩展成 y ~ x + h，h 是一个新变量，\(h=\sum(z+w)^2\) |
function |
将数学函数使用在方程中， log(y) ~ x + z + w，将预测 log(y) 而不是 y |
\[ RSS=\sum^n_{i=1}(y_i - \hat{y_i})^2 \]
fit <- lm(weight ~ height, data=women)
plot(women$height, women$weight)
abline(fit)
summary(fit)##
## Call:
## lm(formula = weight ~ height, data = women)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7333 -1.1333 -0.3833 0.7417 3.1167
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -87.51667 5.93694 -14.74 1.71e-09 ***
## height 3.45000 0.09114 37.85 1.09e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.525 on 13 degrees of freedom
## Multiple R-squared: 0.991, Adjusted R-squared: 0.9903
## F-statistic: 1433 on 1 and 13 DF, p-value: 1.091e-14residuals(fit)## 1 2 3 4 5 6
## 2.41666667 0.96666667 0.51666667 0.06666667 -0.38333333 -0.83333333
## 7 8 9 10 11 12
## -1.28333333 -1.73333333 -1.18333333 -1.63333333 -1.08333333 -0.53333333
## 13 14 15
## 0.01666667 1.56666667 3.11666667fitted(fit)## 1 2 3 4 5 6 7 8
## 112.5833 116.0333 119.4833 122.9333 126.3833 129.8333 133.2833 136.7333
## 9 10 11 12 13 14 15
## 140.1833 143.6333 147.0833 150.5333 153.9833 157.4333 160.8833plot(fit)
anova(fit)## Analysis of Variance Table
##
## Response: weight
## Df Sum Sq Mean Sq F value Pr(>F)
## height 1 3332.7 3332.7 1433 1.091e-14 ***
## Residuals 13 30.2 2.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1vcov(fit)## (Intercept) height
## (Intercept) 35.247305 -0.539880952
## height -0.539881 0.008305861回归系数为3.45显著不同于0（p < 0.001）
多项式回归
fit2 <- lm(weight ~ height + I(height^2), data=women)
summary(fit2)##
## Call:
## lm(formula = weight ~ height + I(height^2), data = women)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.50941 -0.29611 -0.00941 0.28615 0.59706
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 261.87818 25.19677 10.393 2.36e-07 ***
## height -7.34832 0.77769 -9.449 6.58e-07 ***
## I(height^2) 0.08306 0.00598 13.891 9.32e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3841 on 12 degrees of freedom
## Multiple R-squared: 0.9995, Adjusted R-squared: 0.9994
## F-statistic: 1.139e+04 on 2 and 12 DF, p-value: < 2.2e-16plot(women$height, women$weight)
lines(women$height,fitted(fit2))
library(car)
scatterplot(weight ~ height, data=women, spread=FALSE, lty.smooth=2, pch=19, main="Women Age 30-39", xlab="Height (inches)", ylab="Weight (lbs.)")
多元线性回归
states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")])
cor(states)## Murder Population Illiteracy Income Frost
## Murder 1.0000000 0.3436428 0.7029752 -0.2300776 -0.5388834
## Population 0.3436428 1.0000000 0.1076224 0.2082276 -0.3321525
## Illiteracy 0.7029752 0.1076224 1.0000000 -0.4370752 -0.6719470
## Income -0.2300776 0.2082276 -0.4370752 1.0000000 0.2262822
## Frost -0.5388834 -0.3321525 -0.6719470 0.2262822 1.0000000scatterplotMatrix(states, spread=FALSE, main="Scatter Plot Matrix")
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
summary(fit)##
## Call:
## lm(formula = Murder ~ Population + Illiteracy + Income + Frost,
## data = states)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7960 -1.6495 -0.0811 1.4815 7.6210
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.235e+00 3.866e+00 0.319 0.7510
## Population 2.237e-04 9.052e-05 2.471 0.0173 *
## Illiteracy 4.143e+00 8.744e-01 4.738 2.19e-05 ***
## Income 6.442e-05 6.837e-04 0.094 0.9253
## Frost 5.813e-04 1.005e-02 0.058 0.9541
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.535 on 45 degrees of freedom
## Multiple R-squared: 0.567, Adjusted R-squared: 0.5285
## F-statistic: 14.73 on 4 and 45 DF, p-value: 9.133e-08带有相互作用的多元线性回归
library(effects)
fit <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)
summary(fit)##
## Call:
## lm(formula = mpg ~ hp + wt + hp:wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0632 -1.6491 -0.7362 1.4211 4.5513
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 49.80842 3.60516 13.816 5.01e-14 ***
## hp -0.12010 0.02470 -4.863 4.04e-05 ***
## wt -8.21662 1.26971 -6.471 5.20e-07 ***
## hp:wt 0.02785 0.00742 3.753 0.000811 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.153 on 28 degrees of freedom
## Multiple R-squared: 0.8848, Adjusted R-squared: 0.8724
## F-statistic: 71.66 on 3 and 28 DF, p-value: 2.981e-13plot(effect("hp:wt", fit, vcov, list(wt=c(2.2,3.2,4.2))), multiline=TRUE)
回归诊断
上述的 lm 方法，只会根据我们给的公式进行拟合，并产生模型，并没有输出说模型是不是合适。有时还会发生，明明没有关系的输入与响应，最后的结论却是二者有关。
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
confint(fit)## 2.5 % 97.5 %
## (Intercept) -6.552191e+00 9.0213182149
## Population 4.136397e-05 0.0004059867
## Illiteracy 2.381799e+00 5.9038743192
## Income -1.312611e-03 0.0014414600
## Frost -1.966781e-02 0.0208304170以上结果表明，文盲率1%的改变，谋杀率有95%的可能在 \([2.38, 5.90]\) 间。因为 Frost 置信区间包含0，我们可以得出谋杀率和气温没有关系。这些推断只是在数据符合模型下的统计假设时才的出的。
回归诊断的典型方法
fit <- lm(weight ~ height, data=women)
par(mfrow=c(2,2))
plot(fit)
fit2 <- lm(weight ~ height + I(height^2), data=women)
par(mfrow=c(2,2))
plot(fit2)
# 去掉异常值
newfit <- lm(weight~ height + I(height^2), data=women[-c(13,15),])
par(mfrow=c(2,2))
plot(newfit)
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
par(mfrow=c(2,2))
plot(fit)
回归诊断的加强方法
library(car)
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
qqPlot(fit, labels=row.names(states), simulate=TRUE, main="Q-Q Plot")
states["Nevada",]## Murder Population Illiteracy Income Frost
## Nevada 11.5 590 0.5 5149 188fitted(fit)["Nevada"]## Nevada
## 3.878958residuals(fit)["Nevada"]## Nevada
## 7.621042rstudent(fit)["Nevada"]## Nevada
## 3.542929从上我们看出，Nevada 有很高的谋杀率（11.5），但是预测出的值确很小（3.878958）。问题出在哪里？
residplot <- function(fit, nbreaks = 10) {
z <- rstudent(fit)
hist(z, breaks = nbreaks, freq = FALSE, xlab="Studentized Residual",
main="Distribution of Errors")
rug(jitter(z), col = "brown")
curve(dnorm(x, mean = mean(z), sd = sd(z)), add = TRUE,
col = "blue", lwd = 2)
lines(density(z)$x, density(z)$y, col = "red", lwd = 2, lty = 2)
legend("topright",
legend = c( "Normal Curve", "Kernel Density Curve"),
lty=1:2, col=c("blue","red"), cex=.7)
}
residplot(fit)
rstudent(fit)["Nevada"]## Nevada
## 3.542929误差应当符合正态分布。
可以注意到右下角出现了异常值。
# 判断是否有自相关性
durbinWatsonTest(fit)## lag Autocorrelation D-W Statistic p-value
## 1 -0.2006929 2.317691 0.248
## Alternative hypothesis: rho != 0# lag为1，表明每个观测值都是和下一个值进行比较的。
# 判断是否有线性关系
crPlots(fit)
# 同方差性
ncvTest(fit)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 1.746514 Df = 1 p = 0.1863156spreadLevelPlot(fit)
##
## Suggested power transformation: 1.209626推荐的指数 p（\(Y^p\)）
线性模型假设的全局检验
library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)##
## Call:
## lm(formula = Murder ~ Population + Illiteracy + Income + Frost,
## data = states)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7960 -1.6495 -0.0811 1.4815 7.6210
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.235e+00 3.866e+00 0.319 0.7510
## Population 2.237e-04 9.052e-05 2.471 0.0173 *
## Illiteracy 4.143e+00 8.744e-01 4.738 2.19e-05 ***
## Income 6.442e-05 6.837e-04 0.094 0.9253
## Frost 5.813e-04 1.005e-02 0.058 0.9541
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.535 on 45 degrees of freedom
## Multiple R-squared: 0.567, Adjusted R-squared: 0.5285
## F-statistic: 14.73 on 4 and 45 DF, p-value: 9.133e-08
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = fit)
##
## Value p-value Decision
## Global Stat 2.7728 0.5965 Assumptions acceptable.
## Skewness 1.5374 0.2150 Assumptions acceptable.
## Kurtosis 0.6376 0.4246 Assumptions acceptable.
## Link Function 0.1154 0.7341 Assumptions acceptable.
## Heteroscedasticity 0.4824 0.4873 Assumptions acceptable.多重共线性（Multicollinearity）
多重共线性是：多个输入变量组合起来的时候和响应有关系，但是单个输入变量确和响应没有关系。
vif(fit)## Population Illiteracy Income Frost
## 1.245282 2.165848 1.345822 2.082547sqrt(vif(fit)) > 2 # problem?## Population Illiteracy Income Frost
## FALSE FALSE FALSE FALSE以上检测说明我们的输入变量中不存在多重共线性的问题。
不寻常的观测值
异常值
异常值是指那些预测值与实际值偏差很大的值。
outlierTest(fit)## rstudent unadjusted p-value Bonferonni p
## Nevada 3.542929 0.00095088 0.047544高杠杆率点（High leverage points）
高杠杆率点：
high leverage are outliers with regard to the other predictors. In other words, they have an unusual combination of predictor values. The response value isn’t involved in determining leverage.
高杠杆率点可以通过 hat 统计的方法来得到，对于给定的数据集，hat 值的平均值为 \(p/n\) ，其中 p 为模型中估计参数的个数，包括截距，n 为数据集的条目数。大于这个平均值 2 到 3 倍的观测值的 hat 值，都是高杠杆率点。
hat.plot <- function(fit) {
p <- length(coefficients(fit))
n <- length(fitted(fit))
plot(hatvalues(fit), main="Index Plot of Hat Values")
abline(h = c(2,3)*p/n, col = "red", lty = 2)
identify(1:n, hatvalues(fit), names(hatvalues(fit)))
}
hat.plot(fit)
## integer(0)影响点（Influential observations）
影响点是那些对模型参数有不成比例的影响的观测值。例如，去掉某个观测值，模型会发生很大的变动。
Cook’s D 值大于 \(4/(n-k-1)\), k 是输入变量的个数, 就表明是影响点。
cutoff <- 4/(nrow(states)-length(fit$coefficients)-2)
plot(fit, which=4, cook.levels=cutoff)
abline(h=cutoff, lty=2, col="red")
有时候 Cook’s D 为1，比 \(4/(n-k-1)\) 更有用。
Cook’s D 只能标识出影响点，却没有这些点怎样影响模型的信息。 Added-Variable Plots 可以给出这些信息。
avPlots(fit, ask=FALSE, id.method="identify")
coefficients(fit)## (Intercept) Population Illiteracy Income Frost
## 1.2345634112 0.0002236754 4.1428365903 0.0000644247 0.0005813055每幅图中的直线表示该输入变量的系数。
influencePlot 可以画出以上提到的3中点
influencePlot(fit, id.method="identify", main="Influence Plot",
sub="Circle size is proportional to Cook’s distance")
在水平线 +2，-2 上下的点是异常值，垂直线0.2或 0.3右边的点。圆圈的大小表示该观测值对模型影响的大小。
矫正措施
删除观测值
通常是删掉影响点或者异常值，然后再生成模型。这种删除要慎重。
变换变量
当模型违反正态假设的时候，一种是修改响应变量的次数，\(Y\rightarrow Y^\lambda\)。如果 \(Y\) 是一个比值，则可以用 \(ln(Y/1-Y)\)
summary(powerTransform(states$Murder))## bcPower Transformation to Normality
##
## Est.Power Std.Err. Wald Lower Bound Wald Upper Bound
## states$Murder 0.6055 0.2639 0.0884 1.1227
##
## Likelihood ratio tests about transformation parameters
## LRT df pval
## LR test, lambda = (0) 5.665991 1 0.01729694
## LR test, lambda = (1) 2.122763 1 0.14512456从上可以看出，推荐的参数为0.6（\(Murder^{0.6}\)），因为0.6接近与0.5，则可以考虑使用开平方。但是在这个例子中 \(\lambda=1\) 不能被拒绝掉（\(p=0.145\)），则没有强有力的证据，证明需要这样一个变换。
另一种是对输入变量进行变换，
boxTidwell(Murder~Population+Illiteracy,data=states)## Score Statistic p-value MLE of lambda
## Population -0.3228003 0.7468465 0.8693882
## Illiteracy 0.6193814 0.5356651 1.3581188
##
## iterations = 19以上结果表明 \(Population^{0.87}\,,Population^{1.36}\)可以提升线性关系。但二者的 p 值表明没必要进行这种变换。
添加或删除变量
要处理多重共线性，删除变量是一个很重要的方法。当我们只是进行预测，多重共线性对我们没有什么影响。而当我们要理解预测值与每个输入之间对关系时，就需要处理多重共线性。一种方法时删除 \(\sqrt{vif}>2\) 对变量。另一种是使用脊回归。
使用另外的回归方法
选择最好的回归模型
- 在学习模型的时候，是不是该包含所有的变量？
- 应该去掉对预测没有贡献对变量吗？
- 为了提高拟合，应该添加高阶或者相互作用吗？
选择精准的模型，选择简单的模型。
比较模型
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
fit2 <- lm(Murder ~ Population + Illiteracy, data=states)
anova(fit1)## Analysis of Variance Table
##
## Response: Murder
## Df Sum Sq Mean Sq F value Pr(>F)
## Population 1 78.854 78.854 12.2713 0.001052 **
## Illiteracy 1 299.646 299.646 46.6307 1.83e-08 ***
## Income 1 0.057 0.057 0.0089 0.925368
## Frost 1 0.021 0.021 0.0033 0.954148
## Residuals 45 289.167 6.426
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(fit2)## Analysis of Variance Table
##
## Response: Murder
## Df Sum Sq Mean Sq F value Pr(>F)
## Population 1 78.854 78.854 12.813 0.0008122 ***
## Illiteracy 1 299.646 299.646 48.690 8.832e-09 ***
## Residuals 47 289.246 6.154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(fit2, fit1)## Analysis of Variance Table
##
## Model 1: Murder ~ Population + Illiteracy
## Model 2: Murder ~ Population + Illiteracy + Income + Frost
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 47 289.25
## 2 45 289.17 2 0.078505 0.0061 0.9939模型1包含在模型2中，因为非显著性测试为0.9939，则说明二者相同，所以我们第2个模型，即去掉 Income 和 Frost。
Akaike Information Criterion (AIC) 提供了另一种比较模型的方法（考虑进了达到拟合的参数个数）。更小的 AIC 值，意味着较少的参数就达到了较高的精准度，也更适合。
AIC(fit1,fit2)## df AIC
## fit1 6 241.6429
## fit2 4 237.6565需要注意的是 anova 需要嵌套的模型，而 AIC 则不需要。
变量的选择
逐步回归（STEPWISE REGRESSION）
添加／删除变量，直到达到某个标准。例如前向逐步回归是每次增加一个变量，直到模型不再提升；后向逐步回归则是一开是包含所有变量，每次去掉一个变量，直到模型开始降低。逐步回归则是前向+后向。
library(MASS)
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
stepAIC(fit, direction="backward")## Start: AIC=97.75
## Murder ~ Population + Illiteracy + Income + Frost
##
## Df Sum of Sq RSS AIC
## - Frost 1 0.021 289.19 95.753
## - Income 1 0.057 289.22 95.759
## <none> 289.17 97.749
## - Population 1 39.238 328.41 102.111
## - Illiteracy 1 144.264 433.43 115.986
##
## Step: AIC=95.75
## Murder ~ Population + Illiteracy + Income
##
## Df Sum of Sq RSS AIC
## - Income 1 0.057 289.25 93.763
## <none> 289.19 95.753
## - Population 1 43.658 332.85 100.783
## - Illiteracy 1 236.196 525.38 123.605
##
## Step: AIC=93.76
## Murder ~ Population + Illiteracy
##
## Df Sum of Sq RSS AIC
## <none> 289.25 93.763
## - Population 1 48.517 337.76 99.516
## - Illiteracy 1 299.646 588.89 127.311##
## Call:
## lm(formula = Murder ~ Population + Illiteracy, data = states)
##
## Coefficients:
## (Intercept) Population Illiteracy
## 1.6515497 0.0002242 4.0807366<none> 的 AIC 值表示没有变量移除时，模型的 AIC。
第一步中，Frost 被移除，AIC 从97.75下降到95.75。该方法有个限制就是并不会评估所有的模型。而下面要介绍的这个方法讲会克服这个限制。
全集回归（ALL SUBSETS REGRESSION）
library(leaps)
leaps <-regsubsets(Murder ~ Population + Illiteracy + Income + Frost, data=states, nbest=4)
plot(leaps, scale="adjr2")
上图显示，最好的模型是两个变量的模型（Population 和 Illiteracy）。
subsets(leaps, statistic="cp", main="Cp Plot for All Subsets Regression", legend=c(3.5, 50))
abline(1,1,lty=2,col="red")
进一步的分析
交叉验证
k-fold 交叉验证：将样本分成 k 个子样本，每个子样本作为测试集，同时除了本身之外的所有子集会构成训练集。
library(bootstrap)
shrinkage <- function(fit, k=10){
theta.fit <- function(x,y){ lsfit(x,y) }
theta.predict <- function(fit,x){ cbind(1,x) %*% fit$coef }
x <- fit$model[,2:ncol(fit$model)]
y <- fit$model[,1]
results <- crossval(x, y, theta.fit, theta.predict, ngroup=k)
r2 <- cor(y, fit$fitted.values)^2
r2cv <- cor(y, results$cv.fit)^2
cat("Original R-square =", r2, "\n")
cat(k, "Fold Cross-Validated R-square =", r2cv, "\n")
cat("Change =", r2-r2cv, "\n")
}
fit <- lm(Murder ~ Population + Income + Illiteracy + Frost, data=states)
shrinkage(fit)## Original R-square = 0.5669502
## 10 Fold Cross-Validated R-square = 0.494881
## Change = 0.07206925fit2 <- lm(Murder~Population+Illiteracy,data=states)
shrinkage(fit2)## Original R-square = 0.5668327
## 10 Fold Cross-Validated R-square = 0.5281635
## Change = 0.03866918相对重要性
哪个变量在预测响应的时候最重要？
一种是看标准回归系数的大小
zstates <- as.data.frame(scale(states))
zfit <- lm(Murder~Population + Income + Illiteracy + Frost, data=zstates)
coef(zfit)## (Intercept) Population Income Illiteracy Frost
## -2.054026e-16 2.705095e-01 1.072372e-02 6.840496e-01 8.185407e-03再有就是 relaimpo 包中的方法。
relweights <- function(fit,...){
R <- cor(fit$model)
nvar <- ncol(R)
rxx <- R[2:nvar, 2:nvar]
rxy <- R[2:nvar, 1]
svd <- eigen(rxx)
evec <- svd$vectors
ev <- svd$values
delta <- diag(sqrt(ev))
lambda <- evec %*% delta %*% t(evec)
lambdasq <- lambda ^ 2
beta <- solve(lambda) %*% rxy
rsquare <- colSums(beta ^ 2)
rawwgt <- lambdasq %*% beta ^ 2
import <- (rawwgt / rsquare) * 100
lbls <- names(fit$model[2:nvar])
rownames(import) <- lbls
colnames(import) <- "Weights"
barplot(t(import),names.arg=lbls,
ylab="% of R-Square",
xlab="Predictor Variables",
main="Relative Importance of Predictor Variables",
sub=paste("R-Square=", round(rsquare, digits=3)),
...)
return(import)
}
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
relweights(fit, col="lightgrey")
## Weights
## Population 14.723401
## Illiteracy 59.000195
## Income 5.488962
## Frost 20.787442